from __future__ import with_statement
import sys

__all__ = ["Counter", "Scanner"]

class Counter:
    def __init__(self):
        Counter.reset(self)

    def __repr__(self):
       return "Num. Sites: %d\n" \
               "Num. Transitions: %d\tPercent sites with Transitions: %.3f\n"\
               "Num. Transversions: %d\tPercent sites with Transversions: %.3f\n"\
               "Num. Mutations: %d\tPercent sites with Mutations: %.3f\n" % \
               ( self.count, \
               self.nTransition, self.percentTransition, \
               self.nTransversion, self.percentTransversion, \
               self.nMutation, self.percentMutation)

    def reset(self):
        self.nTransition = 0
        self.nTransversion = 0
        self.count = 0

    @property 
    def nMutation(self):
        """The sum of transitions and transversions"""
        return self.nTransition + self.nTransversion

    @property
    def percentTransition(self):
        if self.count == 0:
            return 0

        return (float(self.nTransition) / self.count) * 100
    @property
    def percentTransversion(self):
        if self.count == 0:
            return 0

        return (float(self.nTransversion) / self.count) * 100
    @property
    def percentMutation(self):
        if self.count == 0:
            return 0

        return (float(self.nMutation) / self.count) * 100

def isTransition(baseA, baseB):
    return (baseA == 'c' and baseB == 't') or \
            (baseA =='t' and baseB == 'c') or \
            (baseA == 'g' and baseB == 'a') or \
            (baseA == 'a' and baseB == 'g')

def notValid(base):
    return  not (base == 'c' or base == 'a' or base  == 't' or base == 'g')

class Scanner:
    def __init__(self):
        self.cg = Counter()
        self.cng = Counter()
        self.other = Counter()
        self.notPms = Counter()
        self.seqLength = 0        

    def resetCounts(self):
        self.cg.reset()
        self.cng.reset()
        self.other.reset()
        self.notPms.reset()
        self.seqLength = 0
    
    @property
    def nValid(self):
        return self.notPms.count + self.cg.count + self.cng.count + \
                self.other.count

    @property
    def nInvalid(self):
        return self.seqLength - self.nValid

    @property
    def nTransition(self):
        return self.cg.nTransition + self.cng.nTransition + \
                self.other.nTransition + self.notPms.nTransition

    @property
    def nTransversion(self):
        return self.cg.nTransversion + self.cng.nTransversion + \
                self.other.nTransversion + self.notPms.nTransversion

    @property    
    def nMutation(self):
        return self.nTransition + self.nTransversion

    @property
    def nPms(self):
        return self.cg.count + self.cng.count + self.other.count

    @property
    def nTransitionInPms(self):
        return self.cg.nTransition + self.cng.nTransition + \
                self.other.nTransition

    @property
    def nTransversionInPms(self):
        return self.cg.nTransversion + self.cng.nTransversion + \
                self.other.nTransversion

    @property
    def nMutationInPms(self):
        return self.nTransitionInPms + self.nTransversionInPms

    @property
    def percentMutationInPms(self):
        if self.nPms == 0:
            return 0

        return (float(self.nMutationInPms) / self.nPms) * 100

    @property
    def percentTransitionInPms(self):
        if self.nPms == 0:
            return 0            

        return (float(self.nTransitionInPms) / self.nPms) * 100

    @property
    def percentTransversionInPms(self):
        if self.nPms == 0:
            return 0

        return (float(self.nTransversionInPms) / self.nPms) * 100

    def scan(self, seqA, seqB, start=0, end=-1):
        seqA = seqA.strip()
        seqB = seqB.strip()

        length = len(seqA)
        if length != len(seqB):
            raise RuntimeError("sequences must be of same length")

        if end < 0:
            end = length

        if start > end:
            start, end = end, start

        if end > length:
            end = length

        if start < 0:
            start = 0

        self.resetCounts()
        self.seqLength = end - start
      
        def isCg(i):
            if i + 1 < length:
                if seqA[i:i+2] == 'cg' or seqB[i:i+2] == 'cg':
                    return True
            if i - 1 >= 0:
                if seqA[i-1:i+1] == 'cg' or seqB[i-1:i+1] == 'cg':
                    return True
            return False
        def isCng(i):
            if i + 2 < length:
                if (seqA[i] == 'c' and seqA[i+2] == 'g') or \
                        (seqB[i] == 'c' and seqB[i+2] == 'g'):
                    return True
            if i - 2 >= 0:
                if (seqA[i-2] == 'c' and seqA[i] == 'g') or \
                        (seqB[i-2] == 'c' and seqB[i] == 'g'):
                    return True
            return False
        def isCorG(i):
            return seqA[i] == 'c' or seqB[i] == 'c' or \
                    seqA[i] == 'g' or seqB[i] == 'g'

        seqA = seqA.lower()
        seqB = seqB.lower()        
        for i in range(start, end):
            if notValid(seqA[i]) or notValid(seqB[i]):
                continue
            if isCg(i):
                self.cg.count += 1
                if isTransition(seqA[i], seqB[i]):
                    self.cg.nTransition += 1
                elif not seqA[i] == seqB[i]:
                    self.cg.nTransversion += 1
            elif isCng(i):
                self.cng.count += 1  
                if isTransition(seqA[i], seqB[i]):
                    self.cng.nTransition += 1
                elif not seqA[i] == seqB[i]:
                    self.cng.nTransversion += 1
            elif isCorG(i):
                self.other.count += 1
                if isTransition(seqA[i],  seqB[i]):
                    self.other.nTransition += 1
                elif not seqA[i] == seqB[i]:
                    self.other.nTransversion += 1
            else:
                self.notPms.count += 1
                if isTransition(seqA[i], seqB[i]):
                    self.notPms.nTransition += 1
                elif not seqA[i] == seqB[i]:
                    self.notPms.nTransversion += 1


def run():
   if len(sys.argv) == 3:
       with open(sys.argv[1]) as fileA:
           seqA = fileA.read()
       with open(sys.argv[2]) as fileB:
           seqB = fileB.read()
   else:        
       seqA = raw_input("Enter first sequence: ")
       seqB = raw_input("Enter second sequence: ")

   scanner = Scanner()    
   scanner.scan(seqA, seqB)

   print "Summary:"
   print "\tSequence length: %d"% scanner.seqLength
   print "\tNumber of valid bases: %d" % scanner.nValid
   print "\tNumber of PMS: %d" % scanner.nPms
   print "\tTotal number of Mutations: %d" % scanner.nMutation
   print "\tTotal number of Mutations in PMS: %d" % scanner.nMutationInPms
   print "\tPercent PMS with Mutations: %.3f" % scanner.percentMutationInPms
   print "\tTotal number of Transitions: %d" % scanner.nTransition
   print "\tTotal number of Transitions in PMS %d" % scanner.nTransitionInPms
   print "\tPercent PMS with Transitions: %.3f" % scanner.percentTransitionInPms
   print "\tTotal number of Transversions: %d" % scanner.nTransversion
   print "\tTotal number of Transversions in PMS: %d" % scanner.nTransversionInPms
   print "\tPercent PMS with Transversions: %.3f" % scanner.percentTransversionInPms 

   print "\nPer-site Detail\n"
   print "CG Potentially Methylated Sites:\n%s\n" % scanner.cg
   print "CNG Potentially Methylated Sites:\n%s\n" % scanner.cng
   print "Other Potentially Methylated Sites:\n%s\n" % scanner.other

if __name__ == "__main__":
    run()
